######################################
######################################
######################################
###Robustness Models -- All Disease###
######################################
######################################
######################################
dis.dat$Affected_flus_all1000 <- dis.dat$Affected_flus_all/1000
dis.dat$lagAffected_flus_all1000 <- dis.dat$lagAffected_flus_all/1000
###########################
####Affected Individuals###
###########################
####Cattle
iv2.cat.aff.flus <-  ivreg(Affected_flus_all1000~cattle_meat_tonnes1000000+t+
                        lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode)|
                        comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.aff.flus, diagnostics=T)
#Cluster results
iv2.cat.aff.sum.flus <- coeftest(iv2.cat.aff.flus, vcov = vcovCL(iv2.cat.aff.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected_flus_all1000","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected_flus_all1000","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.aff.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.aff.flus <-  ivreg(Affected_flus_all1000~chicken_meat_tonnes1000000+t+
                          lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode)|
                          comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.aff.flus, diagnostics=T)
#Cluster results
iv2.chick.aff.sum.flus <- coeftest(iv2.chick.aff.flus, vcov = vcovCL(iv2.chick.aff.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected_flus_all1000","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected_flus_all1000","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.aff.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.aff.flus <-  ivreg(Affected_flus_all1000~pigs_meat_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode)|
                         comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.aff.flus, diagnostics=T)
#Cluster results
iv2.pigs.aff.sum.flus <- coeftest(iv2.pigs.aff.flus, vcov = vcovCL(iv2.pigs.aff.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected_flus_all1000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected_flus_all1000","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagAffected_flus_all1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.aff.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.aff.sum.flus, iv2.chick.aff.sum.flus, iv2.pigs.aff.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.aff.flus, iv2.chick.aff.flus, iv2.pigs.aff.flus)
#Corrected weak instruments
iv.main.aff.flus <- data.frame(cbind(iv2.cat.aff.sum.stock.flus,iv2.chick.aff.sum.stock.flus,iv2.pigs.aff.sum.stock.flus))
colnames(iv.main.aff.flus) <- c("Cattle (Affected)", "Chicken (Affected)", "Pigs (Affected)")
iv.main.aff.flus


##############################
###Cause (Poultry and Pigs)###
##############################

####Chicken
iv2.chick.cause.flus <-  ivreg(cause_poultry_flus~chicken_meat_tonnes1000000+t+
                                 lagcause_poultry_flus+logpopulation+logGDPpc+factor(ccode)|
                                 comb.np.use_tonnes1000000+t+lagcause_poultry_flus+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.cause.flus, diagnostics=T)
#Cluster results
iv2.chick.cause.sum.flus <- coeftest(iv2.chick.cause.flus, vcov = vcovCL(iv2.chick.cause.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("cause_poultry_flus","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagcause_poultry_flus","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagcause_poultry_flus+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagcause_poultry_flus+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.cause.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


####Pigs
iv2.pigs.cause.flus <-  ivreg(cause_pig_flus~pigs_meat_tonnes1000000+t+lagcause_pig_flus+logpopulation+logGDPpc+factor(ccode)|
                                comb.np.use_tonnes1000000+t+lagcause_pig_flus+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.cause.flus, diagnostics=T)
#Cluster results
iv2.pigs.cause.sum.flus <- coeftest(iv2.pigs.cause.flus, vcov = vcovCL(iv2.pigs.cause.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("cause_pig_flus","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagcause_pig_flus","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagcause_pig_flus+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagcause_pig_flus+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.cause.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

##Export IV models table to Latex
#Estimates
stargazer(iv2.chick.cause.sum.flus, iv2.pigs.cause.sum.flus)
#Obs and diagnostics
stargazer(iv2.chick.cause.flus, iv2.pigs.cause.flus)
#Corrected weak instruments
iv.main.cause.flus <- data.frame(cbind(iv2.chick.cause.sum.stock.flus,iv2.pigs.cause.sum.stock.flus))
colnames(iv.main.cause.flus) <- c("Chicken (Cause)", "Pigs (Cause)")
iv.main.cause.flus


#######################################
#######################################
###Alternative Explanatory Variables###
#######################################
#######################################

###########################
###Meat Production Index###
###########################

####Cattle
iv2.cat.PI.flus <-  ivreg(flus_all~cattle_meat_PI+t+
                       lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                       comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.PI.flus, diagnostics=T)
#Cluster results
iv2.cat.PI.sum.flus <- coeftest(iv2.cat.PI.flus, vcov = vcovCL(iv2.cat.PI.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_PI ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_PI ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.PI.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.PI.flus <-  ivreg(flus_all~chicken_meat_PI+t+
                         lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                         comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.PI.flus, diagnostics=T)
#Cluster results
iv2.chick.PI.sum.flus <- coeftest(iv2.chick.PI.flus, vcov = vcovCL(iv2.chick.PI.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_PI ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_PI ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.PI.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.PI.flus <-  ivreg(flus_all~pigs_meat_PI+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                        comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.PI.flus, diagnostics=T)
#Cluster results
iv2.pigs.PI.sum.flus <- coeftest(iv2.pigs.PI.flus, vcov = vcovCL(iv2.pigs.PI.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_PI ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_PI ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.PI.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.PI.sum.flus, iv2.chick.PI.sum.flus, iv2.pigs.PI.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.PI.flus, iv2.chick.PI.flus, iv2.pigs.PI.flus)
#Corrected weak instruments
iv.main.PI.flus <- data.frame(cbind(iv2.cat.PI.sum.stock.flus,iv2.chick.PI.sum.stock.flus,iv2.pigs.PI.sum.stock.flus))
colnames(iv.main.PI.flus) <- c("Cattle (PI)", "Chicken (PI)", "Pigs (PI)")
iv.main.PI.flus


#########################################
#########################################
###Disaggregate Instrumental Variables###
#########################################
#########################################

##############
###Nitrates###
##############
####Cattle
iv2.cat.nit.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+t+
                        lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                        nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.nit.flus, diagnostics=T)
#Cluster results
iv2.cat.nit.sum.flus <- coeftest(iv2.cat.nit.flus, vcov = vcovCL(iv2.cat.nit.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.nit.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.nit.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+t+
                          lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                          nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.nit.flus, diagnostics=T)
#Cluster results
iv2.chick.nit.sum.flus <- coeftest(iv2.chick.nit.flus, vcov = vcovCL(iv2.chick.nit.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.nit.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.nit.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                         nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.nit.flus, diagnostics=T)
#Cluster results
iv2.pigs.nit.sum.flus <- coeftest(iv2.pigs.nit.flus, vcov = vcovCL(iv2.pigs.nit.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.nit.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.nit.sum.flus, iv2.chick.nit.sum.flus, iv2.pigs.nit.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.nit.flus, iv2.chick.nit.flus, iv2.pigs.nit.flus)
#Corrected weak instruments
iv.main.nit.flus <- data.frame(cbind(iv2.cat.nit.sum.stock.flus,iv2.chick.nit.sum.stock.flus,iv2.pigs.nit.sum.stock.flus))
colnames(iv.main.nit.flus) <- c("Cattle (Nitrate)", "Chicken (Nitrate)", "Pigs (Nitrate)")
iv.main.nit.flus



################
###Phosphates###
################

####Cattle
iv2.cat.phos.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+t+
                         lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                         phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.phos.flus, diagnostics=T)
#Cluster results
iv2.cat.phos.sum.flus <- coeftest(iv2.cat.phos.flus, vcov = vcovCL(iv2.cat.phos.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.phos.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.phos.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+t+
                           lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                           phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.phos.flus, diagnostics=T)
#Cluster results
iv2.chick.phos.sum.flus <- coeftest(iv2.chick.phos.flus, vcov = vcovCL(iv2.chick.phos.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.phos.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.phos.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                          phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.phos.flus, diagnostics=T)
#Cluster results
iv2.pigs.phos.sum.flus <- coeftest(iv2.pigs.phos.flus, vcov = vcovCL(iv2.pigs.phos.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.phos.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.phos.sum.flus, iv2.chick.phos.sum.flus, iv2.pigs.phos.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.phos.flus, iv2.chick.phos.flus, iv2.pigs.phos.flus)
#Corrected weak instruments
iv.main.phos.flus <- data.frame(cbind(iv2.cat.phos.sum.stock.flus,iv2.chick.phos.sum.stock.flus,iv2.pigs.phos.sum.stock.flus))
colnames(iv.main.phos.flus) <- c("Cattle (Phosphates)", "Chicken (Phosphates)", "Pigs (Phosphates)")
iv.main.phos.flus



##############################
##############################
###Controls and Confounders###
##############################
##############################


#########################
####Imports and Export###
#########################

##Create log values
dis.dat$logcattle.meat.imp_1kusd <- log(dis.dat$cattle.meat.imp_1kusd+1)
dis.dat$logchicken.meat.imp_1kusd <- log(dis.dat$chicken.meat.imp_1kusd+1)
dis.dat$logpigs.meat.imp_1kusd <- log(dis.dat$pigs.meat.imp_1kusd+1)
dis.dat$logcattle.meat.exp_1kusd <- log(dis.dat$cattle.meat.exp_1kusd+1)
dis.dat$logchicken.meat.exp_1kusd <- log(dis.dat$chicken.meat.exp_1kusd+1)
dis.dat$logpigs.meat.exp_1kusd <- log(dis.dat$pigs.meat.exp_1kusd+1)


####Cattle
iv2.cat.impexp.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+t+
                           lagflus_all+logpopulation+logGDPpc+factor(ccode)+
                           logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                           logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                           comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
                           logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                           logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.cat.impexp.flus, diagnostics=T)
#Cluster results
iv2.cat.impexp.sum.flus <- coeftest(iv2.cat.impexp.flus, vcov = vcovCL(iv2.cat.impexp.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc", "logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.impexp.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.impexp.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+t+
                             lagflus_all+logpopulation+logGDPpc+factor(ccode)+logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                             logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                             comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
                             logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                             logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.chick.impexp.flus, diagnostics=T)
#Cluster results
iv2.chick.impexp.sum.flus <- coeftest(iv2.chick.impexp.flus, vcov = vcovCL(iv2.chick.impexp.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.impexp.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.impexp.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                            logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                            comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
                            logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                            logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.pigs.impexp.flus, diagnostics=T)
#Cluster results
iv2.pigs.impexp.sum.flus <- coeftest(iv2.pigs.impexp.flus, vcov = vcovCL(iv2.pigs.impexp.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.impexp.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.impexp.sum.flus, iv2.chick.impexp.sum.flus, iv2.pigs.impexp.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.impexp.flus, iv2.chick.impexp.flus, iv2.pigs.impexp.flus)
#Corrected weak instruments
iv.main.impexp.flus <- data.frame(cbind(iv2.cat.impexp.sum.stock.flus,iv2.chick.impexp.sum.stock.flus,iv2.pigs.impexp.sum.stock.flus))
colnames(iv.main.impexp.flus) <- c("Cattle (Imp./Exp.)", "Chicken (Imp./Exp.)", "Pigs (Imp./Exp.)")
iv.main.impexp.flus


############################
###Local Consumption Only###
############################
#Normalized by 1M
dis.dat$cattle.meat.exp.ton1000000 <- (dis.dat$cattle.meat.exp.ton/1000000)
dis.dat$chicken.meat.exp.ton1000000 <- (dis.dat$chicken.meat.exp.ton/1000000)
dis.dat$pigs.meat.exp.ton1000000 <- (dis.dat$pigs.meat.exp.ton/1000000)
#Create local consumption indicators
dis.dat$cattle_meat_tonnes1000000_local <- dis.dat$cattle_meat_tonnes1000000-dis.dat$cattle.meat.exp.ton1000000
dis.dat$chicken_meat_tonnes1000000_local <- dis.dat$chicken_meat_tonnes1000000-dis.dat$chicken.meat.exp.ton1000000
dis.dat$pigs_meat_tonnes1000000_local <- dis.dat$pigs_meat_tonnes1000000-dis.dat$pigs.meat.exp.ton1000000
#Correct negative error
dis.dat$cattle_meat_tonnes1000000_local <- ifelse(dis.dat$cattle_meat_tonnes1000000_local<0,0,dis.dat$cattle_meat_tonnes1000000_local)
summary(dis.dat$cattle_meat_tonnes1000000_local)
summary(dis.dat$chicken_meat_tonnes1000000_local)
summary(dis.dat$pigs_meat_tonnes1000000_local)

####Cattle
iv2.cat.local.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000_local+t+
                               lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                               comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.local.flus, diagnostics=T)
#Cluster results
iv2.cat.local.sum.flus <- coeftest(iv2.cat.local.flus, vcov = vcovCL(iv2.cat.local.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000_local ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.local.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.local.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000_local+t+
                                 lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                                 comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.local.flus, diagnostics=T)
#Cluster results
iv2.chick.local.sum.flus <- coeftest(iv2.chick.local.flus, vcov = vcovCL(iv2.chick.local.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000_local ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.local.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.local.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000_local+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                                comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.local.flus, diagnostics=T)
#Cluster results
iv2.pigs.local.sum.flus <- coeftest(iv2.pigs.local.flus, vcov = vcovCL(iv2.pigs.local.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000_local ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.local.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.local.sum.flus, iv2.chick.local.sum.flus, iv2.pigs.local.sum.flus)
#Obs and diagnostics
stargazer( iv2.cat.local.flus, iv2.chick.local.flus, iv2.pigs.local.flus)
#Corrected weak instruments
iv.main.local.flus <- data.frame(cbind(iv2.cat.local.sum.stock.flus,iv2.chick.local.sum.stock.flus,iv2.pigs.local.sum.stock.flus))
colnames(iv.main.local.flus) <- c("Cattle (Local)", "Chicken (Local)", "Pigs (Local)")
iv.main.local.flus


#####################
####No Country FEs###
#####################
####Cattle
iv2.cat.nfes.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+t+
                         lagflus_all+logpopulation+logGDPpc|
                         comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.cat.nfes.flus, diagnostics=T)
#Cluster results
iv2.cat.nfes.sum.flus <- coeftest(iv2.cat.nfes.flus, vcov = vcovCL(iv2.cat.nfes.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.nfes.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.nfes.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+t+
                           lagflus_all+logpopulation+logGDPpc|
                           comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.chick.nfes.flus, diagnostics=T)
#Cluster results
iv2.chick.nfes.sum.flus <- coeftest(iv2.chick.nfes.flus, vcov = vcovCL(iv2.chick.nfes.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.nfes.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.nfes.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc|
                          comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.pigs.nfes.flus, diagnostics=T)
#Cluster results
iv2.pigs.nfes.sum.flus <- coeftest(iv2.pigs.nfes.flus, vcov = vcovCL(iv2.pigs.nfes.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.nfes.sum.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.nfes.sum.flus, iv2.chick.nfes.sum.flus, iv2.pigs.nfes.sum.flus)
#Obs and diagnostics
stargazer(iv2.cat.nfes.flus, iv2.chick.nfes.flus, iv2.pigs.nfes.flus)
#Corrected weak instruments
iv.main.nfes.flus <- data.frame(cbind(iv2.cat.nfes.sum.stock.flus,iv2.chick.nfes.sum.stock.flus,iv2.pigs.nfes.sum.stock.flus))
colnames(iv.main.nfes.flus) <- c("Cattle (NFEs)", "Chicken (NFEs)", "Pigs (NFEs)")
iv.main.nfes.flus


####################
####Year Only FEs###
####################
####Cattle
iv2.cat.yfe.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+
                        lagflus_all+logpopulation+logGDPpc+factor(year)|
                        comb.np.use_tonnes1000000+lagflus_all+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.cat.yfe.flus, diagnostics=T)
#Cluster results
iv2.cat.sum.yfe.flus <- coeftest(iv2.cat.yfe.flus, vcov = vcovCL(iv2.cat.yfe.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.sum.yfe.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.yfe.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+
                          lagflus_all+logpopulation+logGDPpc+factor(year)|
                          comb.np.use_tonnes1000000+lagflus_all+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.chick.yfe.flus, diagnostics=T)
#Cluster results
iv2.chick.sum.yfe.flus <- coeftest(iv2.chick.yfe.flus, vcov = vcovCL(iv2.chick.yfe.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.sum.yfe.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.yfe.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(year)|
                         comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.pigs.yfe.flus, diagnostics=T)
#Cluster results
iv2.pigs.sum.yfe.flus <- coeftest(iv2.pigs.yfe.flus, vcov = vcovCL(iv2.pigs.yfe.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ lagflus_all+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.sum.yfe.stock.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.sum.yfe.flus, iv2.chick.sum.yfe.flus, iv2.pigs.sum.yfe.flus)
#Obs and diagnostics
stargazer(iv2.cat.yfe.flus, iv2.chick.yfe.flus, iv2.pigs.yfe.flus)
#Corrected weak instruments
iv.main.nfes.yfe.flus <- data.frame(cbind(iv2.cat.sum.yfe.stock.flus,iv2.chick.sum.yfe.stock.flus,iv2.pigs.sum.yfe.stock.flus))
colnames(iv.main.nfes.yfe.flus) <- c("Cattle (YFEs)", "Chicken (YFEs)", "Pigs (YFEs)")
iv.main.nfes.yfe.flus

###########################
####Country and Year FEs###
###########################
####Cattle
iv2.cat.cyf.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+factor(year)+
                        logpopulation+logGDPpc+factor(ccode)|
                        comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.cyf.flus, diagnostics=T)
#Cluster results
iv2.cat.sum.cyf.flus <- coeftest(iv2.cat.cyf.flus, vcov = vcovCL(iv2.cat.cyf.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.sum.stock.cyf.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.cyf.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+factor(year)+
                          logpopulation+logGDPpc+factor(ccode)|
                          comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.cyf.flus, diagnostics=T)
#Cluster results
iv2.chick.sum.cyf.flus <- coeftest(iv2.chick.cyf.flus, vcov = vcovCL(iv2.chick.cyf.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.sum.stock.cyf.flus <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.cyf.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode)|
                         comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.cyf.flus, diagnostics=T)
#Cluster results
iv2.pigs.sum.cyf.flus <- coeftest(iv2.pigs.cyf.flus, vcov = vcovCL(iv2.pigs.cyf.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc","year")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ factor(year)+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.sum.stock.cyf.flus  <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.sum.cyf.flus , iv2.chick.sum.cyf.flus , iv2.pigs.sum.cyf.flus )
#Obs and diagnostics
stargazer(iv2.cat.cyf.flus , iv2.chick.cyf.flus , iv2.pigs.cyf.flus )
#Corrected weak instruments
iv.main.nfes.cyf.flus  <- data.frame(cbind(iv2.cat.sum.stock.cyf.flus ,iv2.chick.sum.stock.cyf.flus ,iv2.pigs.sum.stock.cyf.flus ))
colnames(iv.main.nfes.cyf.flus ) <- c("Cattle (CYFEs)", "Chicken (CYFEs)", "Pigs (CYFEs)")
iv.main.nfes.cyf.flus 

################
################
###GMM Models###
################
################
library(plm)

##Keep only variables of interest
subset.dat1.flus <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","lagflus_all","logpopulation","logGDPpc","ccode", "year")]) 
#Convert data to pgmm formal
gmm.dat.flus <- pdata.frame(subset.dat1.flus, index=c("ccode", "year"))


#set seed
set.seed(50)

###Cattle
gmm.cattle.flus <- pgmm(flus_all~cattle_meat_tonnes1000000+t+
                          lagflus_all+logpopulation+logGDPpc| 
                          lag(flus_all,2:99),
                        data = gmm.dat.flus, 
                        effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.cattle.flus)

###Cattle
gmm.chicken.flus <- pgmm(flus_all~chicken_meat_tonnes1000000+t+
                           lagflus_all+logpopulation+logGDPpc| 
                           lag(flus_all,2:99),
                         data = gmm.dat.flus, 
                         effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.chicken.flus)

###Pigs
gmm.pigs.flus <- pgmm(flus_all~pigs_meat_tonnes1000000+t+
                        lagflus_all+logpopulation+logGDPpc| 
                        lag(flus_all,2:99),
                      data = gmm.dat.flus, 
                      effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.pigs.flus)


####Export to LaTex
stargazer(gmm.cattle.flus, gmm.chicken.flus,gmm.pigs.flus)

